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The position of an invasion front, propagating into an unstable state, fluctuates because of the 
shot noise coming from the discreteness of reacting particles and stochastic character of the reactions 
and diffusion. A recent macroscopic theory [Meerson and Sasorov, Phys. Rev. E 84, 030101(R) 
(2011)] yields the probability of observing, during a long time, an unusually slow front. The theory 
is formulated as an effective Hamiltonian mechanics which operates with the density field and the 
conjugate "momentum" field. Further, the theory assumes that the most probable density field 
history of an unusually slow front represents, up to small corrections, a traveling front solution 
of the Hamilton equations. Here we verify this assumption by solving the Hamilton equations 
numerically for models belonging to the directed percolation universality class. 
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I. INTRODUCTION 



X 



The Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) 
equation [l[, 



d t q = q- q 2 + d 2 x q , 



(1) 



describes invasion of an unstable state, q(x — > oo, t) = 0, 
by a stable state, q[x — » — oo, f) = 1. This equation 
serves as a fundamental model in mathematical genetics 
and population biology [J Q- Closely related equations 
appear in chemical kinetics |3j, extreme value statistics 
[i[ , dynamics of disordered systems [f| and even particle 
physics @. 

Invasion fronts correspond to traveling front solutions 
(TFSs) of Eq. ©: £jr(x,t) = Qo,c(0» where € = x - ct. 
Qo,c(0 solves the ordinary differential equation 



?0,c 



cQ' . c + Qo,c - Ql, c = . 



(2) 



where the primes stand for the derivative with respect to 
the (single) argument. For a sufficiently steep initial con- 
dition, the solution of Eq. ([T]) approaches at long times 
the limiting TFS of Eq. ©, Q ,2(£)> with the velocity 
Co = 2, see Ref. [7| for a comprehensive review. This 
special value of the front velocity is determined by the 
dynamics of the leading edge of the front, where one can 
linearize Eq. ([T]) around q = 0. In other words, the non- 
linear front solution, as described by Eq. flU, is "pulled" 
by its leading edge, hence the term "pulled fronts" p}, of 
which the FKPP equation (JTJ) is the best studied exam- 
ple. 
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As a mean-field equation, Eq. (JTJ does not account for 
the shot noise coming from the discreteness of particles 
and stochastic character of the particle reactions and ran- 
dom walk. Even when this noise is week, it causes the 
front position to strongly deviate from the mean-field the- 
ory prediction. This deviation has a systematic part — the 
front velocity shift - and a fluctuating part. If N ^> 1 is 
the effective number of particles in the front region, the 
shifted front velocity is (clHlj 



= 2- 



ln 2 iV 



6tt 2 InQnAO 
ln 3 iV 



(3) 



In its turn, the front diffusion coefficient scales as In - N 
fid^l^ . These anomalous properties of noisy pulled 
fronts are very different from the properties of noisy 
fronts propagating into metastable states. In the latter 
case the front velocity shift and the front diffusion coef- 
ficient both scale as 1 /N and are therefore much smaller 

EMS- 

The front diffusion coefficient is determined by typical, 
relatively small fluctuations of the front position. What 
is the probability of large fluctuations? For negative fluc- 
tuations (that is, unusually slow fronts) this question 
was recently addressed in Ref. [l7j in the framework of 
macroscopic fluctuation theory (MFT). This theory can 
be derived, using WKB approximation, from the mas- 
ter equation for microscopic lattice models that involve 
on-site reactions among particles and unbiased random 
walk [l8|, [r| . The MFT can be formulated as a classi- 
cal Hamiltonian field theory which involves the particle 
density field q{x, t) and the canonically conjugate "mo- 
mentum" field p(x,t). It is similar in spirit to the math- 
ematically rigorous MFT of Bertini, De Sole, Gabrielli, 
Jona-Lasinio, and Landim [2(j, see also Refs. |2lTl23 |. de- 
veloped for diffusive lattice gases in the absence of on-site 
reactions. It is also analogous to the Martin-Siggia-Rose 
field-theoretical formalism [23 for continuous stochastic 
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systems. 

Let Vic) be the probability density that the fluctuat- 
ing front moves, during a long time T 3> 1, with av- 
erage velocity c that is smaller than c*. Within the 
MFT framework, the (minus) logarithm of this proba- 
bility density is proportional to the mechanical action 
along a certain phase trajectory q(t),p(i) of the Hamilto- 
nian system, see below. For a given average front velocity 
c, q(x, t) describes the optimal (most probable) density 
history of the fluctuating front. Meerson and Sasorov 
17] considered the set of reactions A <=* 2 A and random 
walk. The crucial conjecture they made (see also Ref. 
i.15.3) was that, apart from boundary layers at t = and 
t = T, the optimal trajectory is a traveling front solution: 
q(x,t) — q(x — ct) and p(x,t) — p(x — ct), of the Hamil- 
ton equations for q and p. The traveling front ansatz re- 
duces these equations to ordinary differential equations. 
The (instanton-type) solutions of these equations were 
found analytically in Ref. [17| for some values of c [most 
importantly for c close, but not too close, to c»] and nu- 
merically otherwise. The resulting lnP(c) turns out to be 
proportional to T, strongly non-Gaussian with respect to 
c and rapidly falling with an increase of 8c = c* — c. The 
iV-dependence of the quantity — T" 1 ln'P(c) undergoes a 
major change depending on 5c: 

\nV(c) f exp (y ) ,^«6c « 
T \ const 7Ve"7fc, < (5c < 1. 

(4) 

The first line of Eq. (0| follows from the second one 
when 8c -C 7r 2 /ln 2 7V. Furthermore, the first line coin- 
cides (up to the numerical pre-factor which is apparently 
non-universal) with prediction from the phenomenologi- 
cal theory of Derrida et al. [nj]. The validity of MFT 
demands that 8c > 2ir 2 In" [13]. Needless to say, 
N must be huge to have these pronounced asymptotic 
regimes. 

It was observed in Ref. [13] that the aforementioned 
properties of In "P(c) also hold for all sets of on-site reac- 
tions which belong to the directed percolation universal- 
ity class: under condition that the system is sufficiently 
close to the characteristic (transcritical) bifurcation of 
the mean-field theory. This is the setting we address 
here. We present strong numerical evidence that the TFS 
is indeed the true optimal history of the particle density 
field, and it yields the leading-order asymptotics of In V 
as described by the scaling relation (fT2|) . 



II. MODEL 

Consider three on-site reactions: branching A A 2A, 
coagulation 2 A A A and decay A A 0. These reactions 
constitute a Markov birth-death process with the birth 
rate X(ni) = Aorij and the death rate /Lt(rij) = ^o n i + 
(o~o/2)ni(ni — 1), where rii is the number of particles on 
site i. Let us define B — Xo/po and K = 2Aq/(To, and 



assume that K ^> 1. The deterministic rate equation for 
the on-site dynamics is 

h = fi Q n(B - 1 - Bn/K). 

When B > 1, n = is a repelling point, and n = 
K(l — l/B) is an attracting point. Adding to these 
reactions symmetric and independent random walk of 
the particles between neighboring sites, with rate con- 
stant Do, we obtain a spatial model that describes in- 
vasion of the unstable state n — by the stable state 
n = K(l—1/B). We assume that the system is close to its 
(transcritical) bifuraction at B — 1, by putting B = 1+8, 
where < 8 <C 1. Under this assumption the model be- 
comes universal: all models belonging to the directed 
percolation universality class, see e.g. Ref. (26^ . behave 
in the same way when properly rescaled. We also assume 
that the hopping rate is sufficiently high: Dq 3> po8, al- 
lowing a continuum description in space. Then the mean- 
field theory of this system is described by the FKPP 
equation ([1]), where q — n/(K8), t and x are rescaled 
as follows: fioSt — > t and xjl — > x, I — [D /(^ 8)] is 
the characteristic diffusion length, and the lattice spac- 
ing is set to one. In their turn, fluctuations of the front 
position are encoded in the MFT e quat ions that can be 
derived in WKB approximation [l8l |26||: 

dtq - q-q 2 +2qp + d 2 x q, (5) 
d t p = -p - p 2 + 2qp ~ d 2 p. (6) 

Equations §5§ and ([6]) are Hamiltonian, with the Hamil- 
tonian 

/oo 
dxw, w = qp(p— 3+1) — 9 x qd x p, (7) 
-oo 

If p(x, t) = 0, Eq. (O coincides with the FKPP equation 
(p}, whereas Eq. ([6]) is obeyed identically. This is the 
noiseless solution. Noisy fronts have a non-zero p{x, t) 
which enables them to move with average velocities less 
than c*. 

The boundary conditions for Eqs. ([5]) and ^ in x and 
t are specified as follows [13 • At x — > — oo there is a 
stationary distribution of the particle density, sharply 
peaked at q = 1. Therefore, we demand 

q(— oo,i) = l and p{— oo,t) = 0, (8) 

which corresponds to the deterministic fixed point {q — 
l,p = 0) of the on-site Hamiltonian Ho(q,p) — qp(p — 
g+1). At x — oo we demand 

?(oo,t)=0, (9) 

whereas p{x, t) must be bounded at finite x. The front 
positions at t = and t = T ^> 1 are specified by kink- 
like particle density profiles, situated at the distance X 
apart; X can be positive, negative or zero. The kinks 
interpolate monotonically between q(x = — oo) = 1 and 
q(x = oo) = and decay sufficiently rapidly at x — > oo. 
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These boundary conditions specify the problem com- 
pletely 

Once the Hamilton equations ([5]) and ([6]) are solved, 
we can calculate the mechanical action S along the phase 
space trajectory q{x, t), p{x, t) and evaluate the probabil- 
ity density V(X,T): 

-i In V(X,T) ~ S(X,T) = [ dx [ dt (pd t q - w) 

N J-oo Jo 

OO rT 

dx I dtq(x,t)p 2 (x,t), (10) 

- oo JO 

where N = K8 2 £. Notice that the effective number of 
particles N (which we assume to be much greater than 
1) does not coincide with the characteristic number of 
particles inside the diffusion length: the latter quantity 
is equal to K St, not KS 2 £. 

Importantly, the same MFT equations (JSJ) , (JBJ) and (|T0|) 
can be derived from the following Langevin equation for 
a noisy FKPP front: 



d t q = q-q 2 + %q+ ^ =±n(x,t) , (II) 

where rj(x,t) is a Gaussian white noise with zero mean 
and unit variance. 



III. FLUCTUATING TRAVELING FRONT 
CONJECTURE 

The fluctuating traveling front conjecture [HI, H3] as- 
sumes that, at T — > oo, S(X, T) has the following scaling 
form: S(X,T) = TF{X/T), so that 



lnV(X,T) = NT F(X/T). 



(12) 



Furthermore, the large-deviation function J- (X/T) is 
mostly contributed to by a TFS of Eqs. (jSJ and ((6]): 
q = q(x — ct) and p = p(x — ct), where c = X/T. This 
TFS solves the coupled ordinary differential equations 

q" + cq' + q-q 2 +2qp = 0, (13) 
p" - cp' + p + p 2 - 2qp = . (14) 

subject to the boundary conditions 

q(— oo) = 1, p(— oo) = 0, and q(po) = 0, (15) 

whereas p(£) must be bounded at finite £ [Hj]. Similarly 
to the model studied in Ref. [l7|, there is a symmetry 
relation among q and p profiles of the TFS: 



p(0 = - 0. 



(16) 



where £o only depends on c. In its turn, Eq. (TlU)) reduces 
to 



1 f°° 
— ]nV(c)~F(c) = J d^q(0A0- 



(17) 



It is straightforward to solve Eqs. (|13l) and (fT4|) numeri- 
cally by using a shooting algorithm [l5j, [17]. Instead, we 
focus here on the important regime of 2 — c <C 1 where an 
analytic perturbation theory can be developed [TtJ • At 
2-c< 1 the action is mostly gathered from the leading 
edge of the front, £ 3> 1, where 5 « 1. In the (almost 
deterministic) "left region" one has \p\ -C 1, whereas in 
the (scarcely populated) "right region" one has q <C 1. 
A closed analytic theory is possible because there is a 
joint region where \p\ <C 1 and q <C 1 simultaneously. 
The matching calculations coincide with those used in 
Ref. |17j . and we obtain and p(£) in terms of two 
overlapping asymptotics: 



q(0 



QoM) . 



^0 - £ » 1, 
C» 1, 



(18) 



and Eq. (fT6|) for p(^). Here 



Co =7r/ % /2 _ 7+1 + 0(^2^). (19) 
In its turn, the large deviation function (fiTf becomes 



3e 



e 2C Q^(C)rfC^ 0.0074 e" 



(20) 



where Qo,2(£) is fixed by the demand that Qo,2(C >• 1) — 
A^e - ^, whereas the e - '' term is absent. When Sc = c* — 
c < vr 2 In -2 iV, Eq. yields the first line of Eq. 

How large must be T for the traveling front asymp- 
totic to give a dominant contribution to the action? For 
c = X/T not close to 2, the criterion is simply T > 1. 
As c approaches 2, the criterion becomes more stringent. 
Indeed, here the q and p fronts of the TFS are shifted 
by a large distance £0 — 7t/a/2 — c 3> 1. In the large re- 
gion between the fronts one has q <C 1 and |p| <C 1, and 
Eqs. ([5]) and ([6]) can be linearized: 



d t q = q + d^q, 
dtp = -p-d 2 p. 



(21) 
(22) 



Equation (|2"Tj) can be solved forward in time starting at 
t = 0, whereas Eq. ([22)1 can be solved backward in time 
starting from t = T. The corresponding initial condi- 
tions q(x,0) and p{x,T) are localized, so their spread 
is described by the diffusion and anti-diffusion terms in 
Eqs. (|2"Tj) and (|22|) . respectively. Therefore, the transient 
time r that it takes for the TFS of the complete problem 
to set in can be estimated from the condition that the dif- 
fusion length ypr is comparable with £o- Therefore, for 
the TFS to give a dominant contribution to the action, 
we must demand T ^> r which leads to 



2-c 



(23) 



Now, the macroscopic theory is valid when there are 
many particles at the leading edge of the front, £ ~ 
£0 — 7r/%/2 — c. This demand boils down to c, - c > 
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27r 2 In 3 N, with c* from Eq. (|3j). At the border of the ap- 

,ir 2 /\n 2 N, 



plicability region of our theory we have 2 - 
and criterion (f23|) becomes T In 2 N. Interestingly, 
In 2 N coincides with the typical relaxation time of fluc- 
tuations contributing significantly to the front diffusion 
in the phenomenological theory of Derrida et al. fill ] . 
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FIG. 1: (a) The computed density profiles q(x,U) (the solid 
lines) and momentum profiles p(x,ti) (the dashed lines) at 
time moments 0, T/4, T/2, 3T/4 and T. The parameters are 
X = 200, T = 120 and L = 250. (b) The shifted density 
profiles q(x + Ai,U) at time moments T/4, T/2, 3T/4 and T. 
The shifts A, are determined from the condition q(x + Aj = 
0, ti) — 1/2. The collapse of the shifted profiles supports the 
fluctuating traveling front conjecture. 



IV. NUMERICS 




FIG. 2: (a) The time dependence of the front positions of q 
(the solid line) and p (the dashed line), defined by the re- 
lations q(x,t) = 1/2 and and p(x,t) — —1/2, respectively. 
The parameters are X = 200, T = 120 and L = 250. The 
velocity of the traveling front is c = 1.65 which is close to 
X/T — 1.66 .... (b) A numerically computed derivative of 
the accumulated action, ds/dt, versus time. 



equation © backward in time from t = T to t = 0. Cor- 
respondingly, it demands mixed boundary conditions in 
time: q at t = and p at t = T. This presents an obsta- 
cle, as the boundary conditions that we specified involve 
the knowledge of q at both t = and t = T. Fortunately, 
it turns out that the boundary condition for q(x, T) can 
be replaced by a boundary condition for p(x, T) that has 
a form of a kink that interpolates monotonically between 
p(x = — oo) = and p(x — oo) = — 1. For simplicity we 
specified 



To test the predictions of the fluctuating traveling 
front conjecture, we solved the time-dependent equations 
Eqs. ([5]) and ^ numerically without any a priori as- 
sumption about the character of solution. In particular, 
the kink-like density profiles at t = and t — T were 
quite different from those predicted by the TFS. We used 
an iteration al gori thm, originally suggested by Chernykh 
and Stepanov [29( for evaluating the probability density 
of large negative velocity gradients in the Burgers tur- 
bulence. Different modifications of this algorithm were 
used for the determination of the optimal history of the 
density field in several diffusive lattice gas settings, with 
and without on-site reactions [H, HH, H3] ■ The al- 
gorithm iterates the diffusion-type equation ([5]) forward 
in time from t = to t = T, and the anti-diffusion-type 



q(x,t = 0) = d(-x), p(x,t = T) = -9(X-x), (24) 

where 8(x) is the Heaviside step function. Note that in 
this case the momentum field p{x, t) is related to the 
density field q(x, t) [27[ by 

p(x,t) = -q(X-x,T-t), 0<t<T. (25) 

We implemented the Chernykh-Stepanov algorithm us- 
ing implicit finite differences in iterations of Eqs. (J5J) and 
© in a finite box |x| < L/2, where L > X. The bound- 
ary conditions at x — ±oo were replaced by the same con- 
ditions at x = ±L/2. Because of the finiteness of the box 
we also needed a boundary condition for p(L/2,t). We 
checked that the results in the bulk are insensitive to this 
condition: p(x,t) always approaches —1 at sufficiently 
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large x and develops a boundary layer at x = L/ 2 to ac- 
commodate the specified boundary condition at x = L/2. 
Therefore, we simply put p(L/2,t) = — 1 in further com- 
putations. 

We found that the implicit realization of the algo- 
rithm is beneficial for convergence of the iterations. This 
is because the implicit realization guarantees that < 
q(x,t) < 1 and — 1 < p{x,t) < for all x and t, once 
these two double inequalities hold at t = and t = T, 
respectively. For the continuous version of the equations, 
the latter property can be proven as follows. Let q > 
tend to zero at finite x. Consider a small vicinity of the 
minimum point of q(x) where d x q = and d^q > 0. Here 
one can neglect all the terms in Eq. ([5]) except the diffu- 
sion term and see that dtq > at the minimum point. As 
a result, q stays positive. Similarly, q cannot reach 1 at 
any finite x. Indeed, suppose that q(x,t) < 1 and tends 
to 1. Then, introducing e[x, t) = 1 — q and neglecting 
polynomial terms of order e and higher, we obtain from 
Eq. © d t e = d^e - 2p. Since p < and d%e > (we 
consider a vicinity of maximum of q which corresponds 
to a minimum of e), we have dte > 0. Therefore, e cannot 
reach 0, and q cannot reach 1. In a similar fashion, one 
can prove that p{x, t) always stays on the interval (—1,0) 
while evolving backward in time. 

The data presented in Figures Q] and [2] were computed 
for X = 200, T = 120 and L = 250. Figure UK shows the 
numerically found spatial profiles of q{x, U) and p(x, ti) at 
different times ti. One can see that, beyond the boundary 
layers at small t (for q) and at t close to T (for p), a 
uniformly translating front develops. This observation 
is confirmed in Fig. [T]d which shows a perfect collapse 
of shifted profiles q(x + A,,t,-) into a single profile. The 
shifts Aj were determined from the condition q(x + Aj = 
0, ti) = 1/2. Figure^ depicts the positions of the q- and 
p- fronts (defined as the points where q = 1/2 and p — 
— 1/2) versus time. Both of them are straight lines with 
the same slope. The resulting traveling front velocity 
c = 1.65 is close to the expected value X/T = 1.66 .... 

A sharp signature of a fluctuating traveling front is a 
linear i-dependence of the accumulated action 



s(t,X,T) 



dx I dt' q(x,t')p (x,t'). 



(26) 



This linear time dependence is verified in Fig. ^jp which 
shows the time derivative of s(t) evaluated numerically. 
As one can see, this time derivative is constant up to 
transients at times close to and to T. 

We also observed fluctuating TFSs for other sets of 
parameters. Two additional examples are shown in Figs. 
|3] and HI Again, depicted are the front positions of q 
and p versus time and the the time derivative of the ac- 
cumulated action s. Constant front velocities (c = f.O 
and c = —1.0 in the respective cases) and constancy of 
ds/dt versus time are clearly seen, confirming the fluc- 
tuating traveling front conjecture. We also checked that 
the results are insensitive to the exact form of the q- and 
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FIG. 3: A numerically computed derivative of the accumu- 
lated action, ds/dt, versus time, for X = 50, T = 49 and 
L = 100. Inset: the time dependence of the front positions of 
q (the solid line) and p (the dashed line) . The traveling front 
velocity is c = 1.0 which is close to X/T ~ 1.02. 
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FIG. 4: Same as in Fig. [3] but for X = - 
L = 100. The traveling front velocity is c 
close to X/T ~ -0.98. 



50. 



T = 51 and 
-1.0 which is 



p-kinks at t = and t = T respectively: under condi- 
tion that the kinks are well localized. Importantly, in all 
numerical examples we did observe a kink-like density 
profile at t — T, thus validating a posteriori the replace- 
ment of the kink-like boundary condition for q(x, T) by 
a kink-like boundary condition for p(x,T), see Eq. (l24l) . 

Figure [5] shows the large deviation function J 7 , corre- 
sponding to the plateau region of ds/dt observed in the 
time-dependent solutions, versus c. Also shown is T from 
Eq. (fTT)) . obtained by numerically solving, by a shoot- 
ing method, the traveling front equations (TT5|) and (|14l) . 
(See Refs. [HI, [ItJ for details of the shooting method.) 
One can see very good agreement which again confirms 
the traveling front conjecture. When c approaches 2, 
T is exceedingly small, see Eq. (|2TJ1) . and the accuracy 
of our time-dependent solution becomes insufficient to 
probe this regime. 

It is crucial that T{c) is a convex function, see Fig. [S] 
As a result, at given X/T, a TFS with c = X/T has a 
lesser action (and, therefore, a higher probability) than 
a front that first moves slower and then faster, or vice 
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pulled noisy fronts gives a full support to the fluctuat- 
ing traveling front conjecture of Ref. [l7|]. We always 
observed that, except at t close to and T, the optimal 
path of the system, as described by the MFT equations, 
has the form of a traveling front of q and p. We also 
confirmed the scaling behavior (fT2|) of the probability 
to observe an unusually slow front. These findings put 
the macroscopic theory of unusually slow pulled noisy 
fronts [l?} on a solid foundation. A major challenge is 
to develop a theory of unusually fast pulled fronts that 
would go beyond the successful phenomenological theory 
of Dcrrida et al. [ill • 



FIG. 5: Triangles: the large deviation function J-(c) corre- 
sponding to the plateau region of ds/dt obtained from the 
time-dependent solutions. Circles: J-(c) from Eq. (|17[) , ob- 
tained by numerically solving the traveling front equations 
(Tl3j and 1(13)1 . 
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